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ABSTRACT 



Context. Supernova Remnants (SNRs) are thought to be the primary candidates for the sources of Galactic cosmic rays. According 
to Diffusive Shock Acceleration theory, SNR shocks produce a power-law spectrum with index s = 2, perhaps non-linearly modified 
to harder spectra at high energy. Observations of SNRs often indicate particle spectra that are softer than that and show features 
not expected from classical theory. Known drawbacks of the standard approach are the assumption that SNRs evolve in a uniform 
environment, and that the reverse shock does not accelerate particles. Relaxing those assumptions increases the complexity of the 
problem, because one needs reliable hydrodynamical data for the plasma flow as well as good estimates for the magnetic field at the 
reverse shock. 

Aims. We show that these two factors are especially important when modeling young core-collapse SNRs that evolve in a complicated 
circumstellar medium shaped by the winds of progenitor stars. 

Methods. We use high-resolution numerical simulations for the hydrodynamical evolution of the SNR. Instead of parametrizations of 
the magnetic-field profiles inside the SNR, we follow the advection of frozen-in magnetic field inside the SNR, and thus obtain the 
B-field value at all locations, in particular at the reverse shock. To model cosmic-ray acceleration we solve the cosmic-ray transport 
equation in test-particle approximation. 

Results. We find that the complex plasma-flow profiles of core-collapse SNRs significantly modify the particle spectra. Additionally, 
the reverse shock strongly affects the emission spectra and the surface brightness. 
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1. Introduction 

It is often assumed that cosmic rays (CRs) with energies be- 
low the knee of the CR spectrum (~ 10 16 eV) should be ac- 
celerated in Supernova Remnants (SNRs). If this is correct, 
then particle spectra should follow a power-law with index, 
s — 2, with an exponential cut-off at very high energy. However, 
the w ealth of recent observ ational data on, e.g., RX J0852.0 - 
4622 (lAharonian et al.ll2007l). R CW 86 dAharonian et alj |2009). 
SN 1006 d Acero et all 1201 Ol) , Cas A ([Acciari et all [2010; 
Abd oetalJ l2010h. and Tycho's SNR (Ucciari et alJ boi U 
Giordano et al.l 120121) provide strong evidence that the particle 



spectra are significantly softer. Besides, the spectral shape is 
not a pure power law with exponential cut-off, as for exam- 
ple in Cas A. In addition, multi-wavelength images of SNRs 
puzzle obser vers and force them to introduce multi-zone emis- 
sion models (lArava & CuJ2010H l emoine-Goumard et al . 2012; 
lAtovan & Dermerll2012l) . 

Several groups have introduced modifi cations to classical 
Diffusive Shock Acceleration (PSA) theory dAxford et al|l977t 
iKrvmskiil Il977t iBel Il978t iBlandford & Ostrikerl 1 19781). an d 
its nonlinear counter part (NDSA) dMalkov & O'C Drury||200ll) . 
Malk ov et al.l (I2012S ) claim that ion-neutral collisions in the SNR 
vicinity may steepen the energy spectrum of particles by one 
power with respect to classical DSA, due to an evanescence of 
Alfven waves that permits the escape of particles in a certain 
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momentum range. lBlasi et al.1 d2012l) discuss the spectral soften- 
ing arising from the energy and momentum transfer of neutrals 
returning to the upstream region. The mean free path of neutrals 
is larger than the length-scale of shock modification on account 
of CR streaming, and hence the effective compression ratio is 
< 4 a nd the p a rticle spectrum soft even for high energy par- 
ticles. Capriolil d2Q12b argues that the upstream scattering cen- 
ters propagate with Alfven velocity which in the presence of 
very strong magnetic-field amplification reduces the compres- 
sion r atio and leads to soft particle spectra (see also lLerche et al.l 
2000). According to this study, the acc eleration efficiency satu- 
rates at around 30%. llnoue et al.l (120101) consider multiple weak 
secondary shocks that appear when the primary forward shock 
propagates through a medium filled with small dense cloudlets. 
The CRs are re-accelerated at these low Mach number shocks, 
resulting in a softer spectrum. The largest and strongest sec- 
ondary shock in the system is the reverse shock (RS) that propa- 
gates through the ejecta. High-resolution radio and X-ray obser- 
vations su pport the notion that particle accelerat i on can also oc- 
cur at RS dGotthelf et alJl200lb iRho et alJl2003 iDeLanev etail 
120021; ISasaki et all 120061; iHelder & VinM 120081) . but only re- 
cently has this possibility been considered in theoretical calcu- 
lation s dZirakashvili & Aharonianll2010l : IZirakashvili & Ptiiskinl 
I20T1 rMezhinskv et al.1 120 12allbh . relaxing restrictin g assump- 
tions about the magnetic field in the RS region (e.g. lEllison et al.l 
2005). The RS accelerates particles of the ejecta and thus pro- 
vides a second population of relativistic particles besides that 
produced at the forward shock. The superposition of these two 
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components may lead to spectral modifications in the volume- 
integrated emission from SNRs that are not easily reproduced 
with so-called one-zone/one-shock models, and so multiple 
zones/particle populations need to be employed (lArava & Cui 
I2010t iLemoine-Goumard et al.ll2012t lAtovan & DermerlfcOlilT 
In any case misinterpretations are possible if the forward shock 
(FS) is considered the sole accelerator. 

Another idealization is the assumption that SNRs evolve 
in a uniform environment. Given that stellar mass-loss con- 
siderably modifies the surroundings of stars, as well as 
known inhomogeneities in the interstellar medium, it is clear 
that this cannot hold, and calculations have started to take 
into a ccount the comp l icated nature of the SNR e n viron- 
ment (lEllison & Bvkovl 12011b IZ irakashvili & Ptuskinl 1201 it 
Ellis on et alJl2012t TTelezhins kv et al. 2012bl) . This increases the 
complexity of the problem, because one requires reliable hydro- 
dynamics of the plasma flow as well as knowledge about the 
magnetic field (MF) at the reverse shock. 

In this paper we show that taking into account both points, 
namely the acceleration of particles at the RS and the complex 
hydrodynamics of SNRs, is particularly important for modeling 
young core-collapse SNRs evolving in a circumstellar medium 
shaped by the wind from the progenitor star. The resulting 
volume-integrated spectrum of particles and their consequent ra- 
diation is significantly different from that of planar-shock calcu- 
lations and dependent on the type of SNR. To accomplish our 
study, we perform high-resolution simulations of the hydrody- 
namical evolution of the SNR with initial and environmental 
conditions representative of Type Ic and Type IIP Supernovae 
(SNe). We consider the transport of frozen-in MF by the plasma 
flow inside the SNR to trace its evolution inside the remnant and 
particularly in the RS region, separately for the radial and the 
tangential field. We model cosmic -ray acceleration by solving 
the cosmic -ray transport equation in test-particle approximation. 
Finally, we calculate the resulting emission from the SNR and 
construct surface-brightness maps in various energy bands. Thus 
we trace the complex particle distribution resulting from particle 
acceleration at both shocks. 



2. Hydrodynamics 

Our goal is to study how the evolution of SNR shock waves into 
the medium created by their progenitor stars modifies the spectra 
and high-energy emission from these objects. Given the impos- 
sibility of exploring the total diversity of paths that can lead to 
a Supernova explosion, we have instead chosen to sample two 
distinctive types of SNe. These illustrate the detailed structure 
of the circumstellar medium into which SNe advance, and also 
how the various features in the spectra, and the very high energy 
emission, relate to the evolutionary properties of the SN. 

We investigate here Type Ic and Type IIP SNe. We have 
chosen these because they show two different regimes of evo- 
lution, and because many SNe, even those of other types, will 
fall somewhere close to or in between these types. For e.g., Cas 
A is classified as a type lib SN, based on spectra obtained from 
light echoes. However, its evo lution in a Red Supergiant (RSG) 
wind (tChevalier & O ishi 2003), within which it is probably cur- 
rently expanding, would be quite similar to the evolution in a 
wind medium for a Type IIP. Thus we hope to give a flavor of 
the spectra and evolution of many different types of young SNRs 
via these calculations. Herein we describe the assumptions that 
went into the simulations, and the evolution of the remnant in 
each case. 
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Fig. 1. Profiles of the density of circumstellar gas before the ex- 
plosion of Type-Ic and Type-IIP SNe. 



In order to simulate the SN expansion into the ambient 
medium, we first need to understand the nature of the ambient 
medium. Since the medium around core-collapse SNe is shaped 
mainly by mass-loss from the progenitor star, we need to take 
into account the evolution of the progenitors of these SNe. 

Type IIP SNe arise from the explosion of RSGs. In gen- 
eral, these stars start as O or B-type stars while on the main 
sequence, with initial masses between about 8 and 30 so- 
lar masses. In the main-sequence stage these stars lose mass 
via radiatively driven winds, with high wind velocities (> 
1000 kms~') with a mass-loss rate of order 10~ 7 M Q yr _1 . 
The interaction of the main-sequence wind with the surround- 
ing (constant-d ensity) interstellar medium gives rise to a wind- 
blown bubble (Weaver et al. 1977; Garcia-Segura et al. I [19961: 
I van Marie et all 120051: lDwarkadasll2005l l2007t iToala & Arthur! 
1201 lUDwarkadas & Rosenbergl2012l) . The standard structure of 
the wind bubble, going outwards in radius from the star, consists 
of a freely expanding wind ending in a wind termination shock, a 
hot, low-density shocked wind medium, a contact discontinuity, 
outer shock and the external medium. If the wind has constant 
parameters, the freely expanding wind will have a density pro- 
file that decreases as r -2 , while the shocked wind will have a 
more or less constant density. It is possible that, depending on 
the surface temperature and number of ionizing photons emitted 
from the star, there is a dense ioni zed (HII) region in side of the 
constant discontinuity (see for e.g.lDwarkadas 201 1). Since this 
does not happen in all cases, we have not taken it into account. 

As the star moves off the main sequence into the RSG stage, 
it grows considerably in size, the wind mass-loss rate increases 
to about 5xl0 _5 M o yr _I while the velocity drops to a low value 



of order 10 km s . This results in a new pressure equilibrium 
being established. The high density (oc M/ y„ ) of the RSG wind 
leads to the formation of a wind region with density almost 4 
orders of magnitude above that of the main-sequence wind. The 
RSG wind can also end in a shell, but we have chosen for sim- 
plicity not to incorporate it. This leads to the density structure 
shown in Fig Q] with an initial high-density wind with density 
decreasing as r~ 2 , followed by a steep drop in density towards the 
main-sequence wind zone, ending in the main-sequence shell. 

If the initial mass of the star exceeds about 30 M (depend- 
ing on which models are used), then the star will not end its life 
as a RSG. It will leave the main sequence to become a RSG, or, 
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Fig. 2. Time evolution of the plasma outflow velocity (top) and density (bottom) as a function of radius for Type-Ic and Type-IIP 
SNRs. 



in stars with initial mass > 50 M , to become a Luminous Blue 
Variable. Following this, it may shed the outer H and perhaps 
He envelope, turning into a Wolf-Rayet (W-R) star. The mass- 
loss rates of W-R stars are somewhat lower than those of RSG 
winds, but their wind velocities are more than 2 orders of mag- 
nitude higher, leading to wind densities that are more than two 
orders of magnitude lower. The high momentum of the winds 
pushes outwards on the RSG shell, breaking it up in t he process 
and mixing the RSG materia l into the W-R wi nd (Dwarkadas 
I2007t iToala & Arthur! |201 U iDwarkadas & Rosenberg! |2012|) . 
This mixed material approaches, and could bounce back from, 
the main-sequence shell. Eventually the system reaches an equi- 
librium situation which in many ways resembles the main se- 
quence bubble (and has almost the same radius, since the main 
sequence shell is very dense and expanding very slowly, on or- 
der 20-50 km s~') but the extent of the freely expanding wind 
and shocked wind, and the radius of the wind termination shock, 
are different. This leads to the situation shown in FigQ] 

In order to explore the situation without having to simulate 
the entire evolution of the progenitor star and its wind medium, 
we have approximated the wind density profile using reasonable 
average parameters. Following this we investigate the evolution 
of the SN shock wave within the medium. We use the VH- 
1 code, a 3-dimensional finite-differenc e hydrodynamic code 
based on the Piecewise Parabolic Method dColella & Woodward! 
Il984l) . Since we want to study the effects of the complicated en- 
vironment, we have assumed similar characteristics for the SN 
density structure in each case, so that differences in the evolu- 
tion result purely from the difference in the ambient medium. In 
each case we assume an ejecta mass of about 5 M and an ex- 
plosion energy of 10 51 ergs. The ejecta density is flat where the 
plasma flow velocity is below a certain value, u //, and decreases 



as a power law with radius, p el ; oc r 9 , where the flow velo city is 
above u fl dChevalier & Fransson|[T994t lDwarkadasll2005l) . The 
interaction of the ejecta with the wind medium sets up a double- 
shock structure as expected, consisting of a forward and reverse 
shock separated by a contact discontinuity. The region between 
the outer shock and contact discontinuity contains shocked sur- 
rounding medium, whereas the region between the inner shock 
and contact discontinuity includes shocked ejecta. The struc- 
ture of this shocked region depends on the ej ecta density and 
the density profile of the surrounding medium dCheva lier 1982; 
lDwarkadasll201 ll) . and has implications for the emission from 
the remnant. 



The shock expansion in the wind medium in both cases is 
quite similar initially, except that the shock moves slower in the 
high-density region, as expected. The evolution changes once 
the shock reaches the end of the freely expanding wind region. 
In the RSG case one finds a huge drop in density beyond the 
wind region, whereas in the W-R case one finds an increase in 
density by a factor of 4. This change in density and the transi- 
tion from the wind zone to a medium of almost constant density 
destroys the self-similarity of the solution. The interaction of the 
SNR forward shock with the wind termination shock leads to a 
reflected shock travelling back into the ejecta in the case of the 
Type Ic SN. In the case of the Type IIP, the steep drop leads to 
the formation of a complicated ejecta structure. These structures 
are used to compute the acceleration of particles at the shock 
fronts. The evolution of the plasma velocity and density profiles 
is shown in Fig. [2] 
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3. Magnetic Field 

The acceleration and subsequent radiation of relativistic par- 
ticles significantly depends on the distribution of MF in the 
SNR, especially in the upstrea m and downstream vicini ties of 
the shocks. In our earlier work dTelezhinskv et al.ll2012alfbl) . we 
assumed parametrizations for the MF inside the remnant follow- 
ing either pressure or density distributions. We found that these 
simple parametrizations fail for Type-IIP SNRs, since the energy 
contained in MF would be roughly equivalent to the explosion 
energy after a few hundreds years of evolution. In this work we 
implement a more realistic description of the MF distribution in- 
side SNRs, which is especially important for the acceleration of 
particles at the RS, a key ingredient of our studies. 

We start with the assumption that the SNR is filled with a 
perfectly conducting fluid, which is common in astrophysical 
plasmas. The evolution of the frozen-in magnetic field satisfies 
the following equation 



dB 

dt 



V x (u x B) 



(1) 



where B is the MF and u is the flow speed. Rewriting the cross 
product in Eq.Q]in terms of dot products gives 



^5 = u(v • B) - B(v • u) + (B • v)u - (u • v)B (2) 

at 

Since we have spherically symmetric flow, u r = u r , ug = 0, 
= and 8u r /d9 = 0, du r /d(f> = 0. Setting 9 = tt/2, i.e. Bg and 
Z?^ scale in the same fashion and together form the tangential 
field, and using 



w r n 1 d(r2B, ) x 1 dBe ^ 1 dB * 

V D — U = -rr 1 1 — 

r l r r 88 r od> 



(3) 



we reduce Eq.[2](see Appendix) to a set of expressions describ- 
ing the evolution of the radial and tangential components of B. 



OB,. 



dr 

-i(Bg Ur )- 



Bg'± 

- Us— 



(4) 



Eq. [4] represents the transport of the MF components with the 
flow. To solve them we need to apply initial values and boundary 
conditions. 



Estimates in the literature for Type-lib SN ( Marti- Vidal et al 
20111) and Type-la SN dChomiuk etaf] l2012t iMargutfi et al 



2012), suggest that the MF in very young remnants is of the or- 
der of 50 G at an age of a few days, when the SNR has a radius 
of roughly /?sn = 10 15 cm. To ensure that initially V ■ B = 
in the ejecta, we assume that B e j 4 is flat where u r < u/i and 
fiej,r(/) oc 1/r 2 where u r > Ufl. Likewise we do so for the tan- 
gential field, and scale all components so that the total volume- 
averaged MF strength is ^ 50 G when /?sn = 10 15 cm. As initial 
condition we thus use 



Bi(r,to) = B e] (r) 



R SN 



flsNR(f()) 2 



(5) 



where to is the initial time of our calculation. 

In the current work we limit ourselves to scenarios with no 
amplification of the MF upstream of the FS, and only the CSM 
field is assumed to be transported through the shock to the down- 
stream region. The magnetic field in the circumstellar medium 
(CSM) of Type-Ic and Type-IIP SNRs is shaped by the winds 



of the progenitor stars. We assume that turbulence equalizes all 
components of the MF in the wind zones until they are trans- 
ported through the FS. The profile of the MF in the environment 
of a core-collapse SNR in our description is then given by 



„, , _ ( B h (r) _R st <r<R h 
(r> \fi b (fl b )VIT R b <r<R sh 



(6) 



where B\,(r) = B st R st /r, B st is the MF at the surface of the pro- 
genitor star (* 100 G for WR stars, * 1 G for RSGs), R st is 
the radius of the progenitor star (« 8/? Q for WR star, m 6QQR Q 
for RSG), Rb is the radius of the wind-blown bubble in our sim- 
ulations (« 7 pc for WR star, w 2 pc for RSG), R s h is the ra- 
dius of the swept-up shell (« 30 pc in both cases), and r is 
the distance from the star. Note that the 1/r scaling of the field 
strength results from the advection of stellar magnetic field in 
a wind, colloquially known as the Parker spiral, where it ap- 
plies to the tangential components only. The radial field in the 
wind zone is suppos ed to arise from turbulence, w hich may in 
fact also amplify it dO'C. Drurv & Downesll20l2l) . and should 
be understood as RMS amplitude. Magnetic-fie ld amplification 
by cosmic-ra y streaming ( Lucek&Bell 2000), upstream dy- 
namo action (Beresnvak et al.ll2.009h . the Richtmeyer-Meshkov 



instability (Sano et al. I2012L or vorticity-generation at shocks 



dGiacalone & J okipii 2007) are not considered here 

The solution of Eq. |4] is obtained numerically, using the 
initial conditions described above and the parameters of the 
SNR plasma flow from the simulations described in Section [2] 
The time-dependent solution for different SNR types is plot- 
ted in Fig. [3] Our soluti ons for the MF evolution are consistent 
with earlier 1-D studies ( Rosenau & Franke nthal 1976), but 2-D 
MHD solutions (ISchure et al.ll2009l) show significant distortions 
at the contact discontinuity due to the growth of the Rayleigh- 
Taylor (or Kruskal-Schwarzschild) instability. 



4. Particle Acceleration Method 

Our method is based on time- dependent kinetic calcu la- 
tions in test-particle approximation dTelezhinskv et al .1120 1 2allbh . 
They are implemented by numerical solution of the diffusion- 
advection equation for the differential particle number density 
on a grid co-moving with the shock wave and in spherically sym- 
metric geometry. To resolve the diffusion length of the lowest- 
energy particles, we perform a coordinate transformation that in- 
creases the spatial resolution near the shock. The spatial coordi- 
nate, x, is related to the new coordinate, x*, for which a uniform 
grid is used, by the equation 



(7) 



where /?sh is the shock radius. Thus, with modest resolution in 
the x*, one may achieve very fine resolution in x where it mat- 
ters, namely in the shock region where the newly injected low- 
momentum particles enter the acceleration process. The other 
benefit is the significant extension of the grid toward x » 1 at 
very low computation cost, since the range in x is proportional 
to the third power of x*. The boundary condition for large x then 
does not really affect our solution, since particles do not leak 
out of the grid but rather are distributed over the huge upstream 
volume according to the diffusion properties of the media. We 
include an exponential transition from Bohm to Galactic diffu- 
sion at around 2/?fs- A large Galactic diffusion coefficient also 
compensates for the deteriorating resolution of the spatial grid 
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Fig. 3. Time evolution of the total, tangential and radial magnetic field as a function of radius for Type-Ic (top) and Type-IIP (bottom) 
SNRs. 



beyond 2R fs- Because our method is based on the test-particle 
approximation, we do not permit the CR pressure at the shock to 
exceed 10% o f the ram pressur e. We use a thermal-leakage in- 
jection model (Blasi et al. 2005), and adjust the injection so that 
the CR-pressure limit is not violated. The injection coefficient, a 
free parameter, is taken to be approximately 5 ■ 10~ 6 for Type-Ic 
and for 5 • 10~ 8 Type-IIP SNR. 



5. Results and discussion 

Based on the hydrodynamic simulations described in Section|2j 
and the magnetic-field profiles obtained in Section [3] we com- 
pute time-dependent particle distributions accelerated by the for- 
ward and reverse shocks. We calculate the resulting emission 
due to synchrotron, inverse Compton and pion-decay processes, 
build corresponding surface-brightness maps and discuss the ob- 
servational implications. We present snapshots at the age of 400, 
1000, and 1975 years to trace spectral and morphological evolu- 
tion of emission coming from Type-Ic and Type-IIP SNRs. 



5.1. Particle spectra 

Particle spectra are given in Fig.|4]for Type-Ic and in Fig.|5]for 
Type-IIP SNR. Each plot shows two sets of spectra. The first 
is based on the assumption that the cosmic -ray scattering cen- 
ters do not significantly move relative to the plasma flow. The 
second set, calculated under the assumption of Alfvenic drift of 
scattering centers, is marked with subscript "A". For each we 
show spectra produced by acceleration at the forward and re- 
verse shocks plotted along with a volume-integrated total spec- 
trum from the whole SNR. 



At an age of 400 years, the Type-Ic SNR is still propagat- 
ing through the freely expanding stellar wind. The flow veloc- 
ity and density profiles are given by the classical self-similar 
solutions. The spectrum of protons accelerated by the FS is 
well described by a power-law with exponential cutoff, in agree- 
ment with classical DSA solutions. The magnetic field at the 
FS is still strong at this age, and so particles reach high en- 
ergies. The spectrum of RS-accelerated protons is also consis- 
tent with standard DSA. Note that the volume integration is per- 
formed over both the upstream and the downstream region of 
the RS, whereas the upstream region of the FS is not covered 
(for d etails we refer the reader to Section 5.1.1 of our earlier 
paper (Telezhi nsky et al.ll2012ah ). The magnetic field at the RS 
is considerably smaller than at the FS, and so particle energies 
reach a TeV only. However, the number of particles injected 
into the acceleration process, and hence the density, is larger. 
Therefore, below a few TeV the CR number density at the RS is 
larger than that at the FS, and the total volume-integrated proton 
spectrum displays a downward step around a few TeV. In phe- 
nomenological studies of SNRs such spectra may be interpreted 
as sum of different popul ations of particles coming fro m differ- 
ent zones of the remnant dArava & Cuill2010trAtovan & Permed 
I2012L iLemoine-Goumard et al .1 120121) . Since the MF strength is 
not large enough at either of the shocks to cause significant syn- 
chrotron losses, electron spectra look similar to those of protons. 

The Type-IIP SNR at the same age also evolves in the freely- 
expanding progenitor wind, though this wind has a much higher 
density. Therefore both electron and proton spectra have a much 
higher intensity compared to the Type-Ic SNR. The density at 
the RS of the Type-IIP SNR is lower than at the FS, and so the 
CR intensity is lower. Nonetheless, the contribution of the RS to 
the overall particle spectrum is noticeable. The magnetic field at 
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Fig. 4. Time evolution of the proton (top) and electron (bottom) spectra for Type-Ic SNR. Total volume-integrated spectra (S) are 
the sum of the contributions from the forward (F) and reverse (R) shocks. The model accounting for Alfvenic drift is marked with 
subscript "A". 
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the RS of the Type-IIP SNR is slightly larger than in the Type-Ic 
case because the radius is smaller, and particles are accelerated 
to somewhat higher energies, cutting off at around a few TeV. 
For the same reason, electron losses in the downstream region of 
the RS are more efficient and steepen the electron spectra beyond 
a characteristic energy, E ra i, downstream of the RS, whereas in 
the upstream region of the RS the MF, and hence the electron 
energy losses, are weak. Since the volume-integrated spectrum 
of particles produced by the RS includes both the downstream 
and the upstream contributions, we observe a bump located at 
£max of the RS upstream electron distributi on, which is lo cated 
beyond E ra( j (see Section 5.1.1 of Telezhins ky et al.l d2012al) for 
details). 

At 1000 years, Type-Ic SNR particle spectra produced by the 
FS differ significantly from those expected via standard DSA. 
The significant spectral modification can be traced back to the 
interaction of the FS with the wind termination shock, at around 
700 years. The effective compression ratio at the FS as seen by 
the CRs increased significantly. The acceleration was boosted 
for a short time, and since the interaction time was brief and the 
spectral index hard, a bump at a few hundred Ge V formed, which 
slowly moved towards higher energies as the particles in it were 
further accelerated by the FS. 

By the time of 1000 years, a break at around 1 TeV is still 
clearly visible. When the FS collided with the wind termina- 
tion shock, a reflected shock that propagates inwards formed. 
By coincidence, the time when the reflected shock reaches the 
RS of the remnant is very close to 1000 years. So the RS spec- 
tra shown are rather specific for this particular time, however we 
show them to illustrate the impact of the hydrodynamics on par- 
ticle spectra. In analogy with the FS, the compression ratio at the 
RS increased for a very brief time of interaction, the spectral in- 
dex became very hard, and a bump at a few tens of GeV formed 
and shifted to a few hundred GeV by the time of 1000 years. 

Compared to the Type-Ic, Type-IIP SNRs do not show such 
a strong variation in particle spectra at the same age. The for- 
ward shock of the Type-IIP SNR encountered the MS wind zone 
at around 600 years. Since the shock entered a very dilute envi- 
ronment, no significant additional compression occurred at the 
FS, but rather a very brief de-compression that did not affect CR 
spectra. After the encounter, another reverse shock formed and 
propagated through the FS-compressed CSM, and so we see a 
three-shock structure: the FS propagating through the MS wind 
zone, the internal RS going through the ejecta, and another re- 
verse shock between them propagating through the dense CSM 
compressed earlier by the FS. Although the divergence of the 
plasma flow velocity is negative between the forward and in- 
termediate shock, thus permitting re-acceleration of particles in 
this region and at the intermediate shock itself, we cannot see 
the effect in volume-integrated spectra. The FS speed increased 
when it entered the dilute wind zone. The maximum energy of 
CRs, E max , increased as expected, but the rate of particle in- 
jection dropped significantly. The separation between the for- 
ward and reverse shocks started to increase on account of the 
boosted velocity of the FS, and so the volume enclosed by the 
RS constitutes a diminishing fraction of the total SNR volume. 
Nevertheless, the contribution of the RS-accelerated particles 
remains visible because the RS propagates through very dense 
ejecta. 

Nearly 1000 years later, at the age of 1975 years, the ve- 
locity profiles of both remnant types have relaxed. The multiple 
weak reflected shocks in Type-Ic SNR have dissipated, but the 
CR spectra still retain their effects, though in general they appear 
smoother than before. The proton distributions from both the FS 



and the RS are softer, and the contribution from the RS to the to- 
tal proton spectra remains visible. The electrons accelerated by 
the FS have rather soft spectra as well, while those accelerated 
by the RS dominate the total electron spectra of Type-Ic SNR. 

In the Type-IIP SNR, at an age of 1975 years, we see an ex- 
cess of high-energy particles accelerated by the FS. The electron 
spectra differ only slightly from the proton spectra due to weak 
synchrotron losses. Once the FS entered the dilute zone, the par- 
ticle injection rate dropped four orders of magnitude. It takes 
time for the particle spectra to adjust to the new equilibrium, 
particularly so at very high energies. For the high-energy parti- 
cles the characteristic diffusion-advection timescale, = k/v 2 s , 
is long, and so they continue to diffuse around the shock and gain 
energy while the lower-energy particles escaped from the shock. 
We can consider the bump as a signature of the delayed reaction 
of the high-energy particles to a sudden change in the injection 
rate. 

The impact of Alfvenic drift of the scattering centers is more 
pronounced in Type-Ic SNR spectra than in Type-IIP. Although 
the magnetic field at the FS is roughly equal for both SNR 
types, and at the RS is stronger in Type-IIP SNR, the large 
density at both shocks of Type-IIP renders the Alfven velocity 
low, and so its impact on particle spectra is hardly visible. In 
the SNR of Type Ic, on the other hand, the spectra are notice- 
ably softer than ones calculated without Alfvenic drift taken into 
account. Alfvenic Mach number is rather moderate for Type-Ic 
SNR shocks. 

5.2. Radial distributions. 

The radial distributions of 2-TeV and 20-TeV particles are pre- 
sented at Fig. [6] There is a clear-cut difference between particle 
profiles depending on both energy and SNR type. We show pro- 
files for the models without Alfvenic drift of scattering centers. 
The impact of a finite Alfven velocity on the radial profiles of 
CRs is marginal. Essentially one only observes a slightly lower 
CR intensity in Type-Ic SNRs. 

The 20-TeV particle profiles comprise only particles ac- 
cedlerated at the FS, because the RS does not accelerate CRs 
up this energy on account of the small MF. On the other hand, 
the RS region is prominent in the 2-TeV particle profiles be- 
cause at this energy CRs do not efficiently propagate away from 
their acceleration sites. Since cosmic-ray electrons lose energy 
much faster than protons, their intensity falls very quickly with 
increasing distance from the shocks. 

At the age of 400 years, Type-Ic SNRs show an increased 
number density of 20-TeV particles between the two shocks. 
There, the MF is higher than in other parts of the SNR, and 
so particles are retained. In the ejecta region, where the MF is 
low, particles are distributed rather uniformly with an intensity 
approximately an order lower than in the shocked region. The in- 
tensity of 2-TeV particles shows an enhancement in the region of 
the RS on account of the RS acceleration, but only marginally so 
because the cutoff energy of particles accelerated at the RS is be- 
low 2 TeV. A steep decrease of intensity toward the SNR center 
arises because the diffusion time of 2-TeV particles is large and 
the age of the SNR is insufficient to distribute 2-TeV particles 
uniformly. 

The profiles of 20-TeV particles in Type-IIP SNRs of the 
same age behave similarly. The intensity of 2-TeV CRs shows 
a significant peak at the RS, though, because the MF at the re- 
verse shock of Type-IIP SNR is larger than in Type-Ic SNR. The 
maximum energy of RS-accelerated particles is therefore con- 
siderably higher than in the case of Type-Ic SNRs, in fact it is 
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Fig. 6. Time evolution of the radial distribution of protons (top) and electrons (bottom) at energy of 2 TeV and 20 TeV for different 
SNR types and models. The 2-TeV distributions were downscaled by a factor 100 for better visualization. 



well beyond 2 TeV, thus leading to a prominent contribution of 
CRs from the RS. 

As has been already mentioned, the age of 1000 years hap- 
pens to be a special phase in the evolution of Type-Ic SNR for 
our choice of parameters, and it is interesting to see what hap- 
pens to particle distribution just after the shocks collide. The 
profiles of 20-TeV particle are largely unaffected by the recent 
interaction. The MF at the RS is insufficient to re-accelerate par- 
ticles up to 20 TeV and raise the intensity at this energy. A small 
peak in between two shocks can be explained by the increased 
MF at the contact discontinuity behind one of the weak reflected 
shocks; a similar structure is visible in the profiles of 2-TeV par- 
ticles. This MF creates a barrier for particle diffusion to the SNR 
center. The 2-TeV CR profiles are strongly affected by the recent 
shock collision. An excess of 2-TeV particles is clearly seen at 
the RS and, in fact, also in the volume-integrated particle spectra. 
Since the particles have no time to diffuse around the accelera- 
tion region, the peak is very sharp. 

After 1000 years, the radial distribution of CRs in Type-IIP 
SNR has changed dramatically compared to a few hundred years 
earlier when the FS entered the MS wind zone. At the age of 
1000 years the bulk of the CRs is relatively old, having been ac- 
celerated when the FS was propagating through the dense RSG 
wind. The high intensity of CRs, which were advected inside the 
SNR, joins smoothly with the low intensity of particles freshly 
accelerated by the FS. The profiles of 20-TeV particles are sig- 
nificantly smoothed on account of fast diffusion, whereas the in- 
tensity of 2-TeV particles shows an increase towards the RS. The 
RS propagating through dense ejecta at 0.4/?fs makes a strong 
contribution to the profile of 2-TeV particles. 

After 1975 years, the particle profiles in Type-Ic SNR have 
lost memory of the reflected shocks present earlier. The spa- 
tial distribution is rather smooth with enhanced intensity in the 
shocked region, at 2 TeV particularly so at the RS suggesting 



that acceleration is still ongoing. Low-energy particles are far 
more abundant in the SNR interior than they were at early times. 
At 20 TeV the intensity is nearly constant in the SNR interior. 

The radial profiles of particles in Type-IIP SNR have not 
changed much by from 1000 to 1975 years. The profiles be- 
came smoother, but the peak at the RS at around 0.3/?fs is 
still observed for 2-TeV particles. Interesting is that the in- 
termediate shock (see previous subsection) is now visible as 
an increase in low-energy particle intensity at around 0.9/?fs- 
Although no impact of the intermediate shock is seen in volume- 
integrated spectra, it clearly affects radial profiles, suggesting 
some re-acceleration occurs between the forward and interme- 
diate shocks in addition to particle trapping at the contact dis- 
continuity. This effect is not so relevant for 20-TeV particles 
on account of fast diffusion, and so their distribution is rather 
smooth. 



5.3. Nonthermal emission 

We consider three radiation processes of non-thermal particles: 
Synchrotron and inverse Compton (IC) emission of electrons, 
and neutral pion decays originating in collisions of CR protons 
with protons at rest. Only primary electrons are considered for 
the leptonic processes. The calculated photon distributions are 
plotted at Fig.|7]for Type-Ic SNR and at Fig.[8]for Type-IIP SNR. 
We also calculate, and discuss in context with the radiation spec- 
tra, intensity maps for the emission at characteristic radio, X-ray 
and gamma-ray wavelengths. The maps are shown in Fig. [9] We 
note that in all calculations of volume-integrated spectra or in- 
tensity maps the radial distributions of CRs, the target material 
and magnetic field are taken into account. The result is signif- 
icantly different from simple toy-model calculations assuming 
uniform particle, density and MF profiles inside the SNR. We 
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first give brief details on the method of calculation and continue 
with a discussion of results. 

The relevant MF component for synchrotron emission is that 
in the plane-of-sky, B ± , which is perpendicular to the line of 
sight (LOS). If all components of the MF are equal, then it is 
straightforward to show that the effective value of B ± 0.8B. 
When properly accounting for the transport of the MF to the 
downstream region of the FS, it is clear that the radial and the 
tangential components are not equal. We therefore derive the 
plane-of-sky component of the MF as 

B ± - - s jB^sin 2 ^cos 2 a + B 2 cos 2 % + B 2 sin 2 a + B t B r sin£sin(2a) 

(8) 



where B r is the radial component of the MF, B, = ^jB 2 + B 2 

is the tangential component of the MF, a is the angle between 
the radial direction and the LOS, and £ is the angle between the 
direction of the tangential field and the z-axis. Since £ is un- 
known and in general may assume any value in the range [0,27r], 
there are two options to choose it. One option is to keep it con- 
stant at £ = 7r/4, the other one is to randomly vary it. The latter 
case would assume stro ng turbulence and little radio polarization 
(IStroman & Po hl 2009). Here we choose £, — n/4 and apply the 
sta ndard formula for sy nchrotron emission as for instance used 
bv lStorneretall (119971) . 

The IC em ission is calculated using th e full Klein-Nishina 
cross section (Blumenth al & Gouldlll970l) for relativistic elec- 
trons following ISturner et al. I (11997b . As target photon field for 
the IC scattering we consider the microwave background only. 
We calculate hadronic gamma-ray emissi on according to the 
procedure described by Huan g et al.l d2007l) . 

Thermal emission is not considered here. However, we have 
made its estimates given the SNR parameters under consid- 
eration. Thermal emission should be very bright for young 
Type-IIP SNR on account of the high gas density in the RSG 
wind. However, at very early times Type IIP SNe tend to 
show the lowest X-ray inten sity of all SN types (see Fig. 3 of 
iDwarkadas & Gruszkoll2012l) . The reason may be that the high 
density also provides strong absorption. 

Most of the 1-TeV IC and pion-decay emission is created by 
particles of roughly 20-TeV energy whereas the energy of elec- 
trons producing 1.4-GHz and 3-keV synchrotron maps depends 
on the MF, which is varying with age, location within the SNR, 
and type of SNR. The general trend to be noted from Fig. [7] 
is that the high-energy emission from Type-Ic SNRs is domi- 
nated by IC emission, with a non-negligible contribution from 
pion-decay only at the early age of 400 years. Fig. |8]shows that 
the high-energy emission of Type-IIP SNR is mostly hadronic, 
with a non-negligible contribution from IC only at the late age 
of 1975 years. This is not surprising, because the very thick 
RSG wind provides a good target for hadronic CR interactions in 
Type-IIP SNR, whereas Type-Ic SNRs expand in the low-density 
wind of the W-R star. While the high-energy flux from Type- 
Ic SNRs is increasing with time, that from Type-IIP SNRs is 
decreasing. Type-IIP SNRs are always brighter in gamma-rays 
than Type-Ic remnants with the exception of late times when the 
fluxes are of the same order. The synchrotron emission of both 
SNR types shows the opposite trend. To be noted from Fig.|7}{9] 
is that the spectra and intensity maps have a few distinct parti- 
cle populations and emitting zones. Thus complex spectra and 
morphology arise naturally, but evolve in time. 

At the age of 400 years, the hadronic emission from Type-Ic 
SNRs comes from two regions: the region of dense ejecta pro- 
vides a high-energy bump in the spectrum, and the region around 



the contact discontinuity accounts for a low-energy bump. The 
energetic CRs accelerated at the FS penetrate deep into the SNR 
and illuminate the ejecta, while low-energy CRs are trapped 
around the RS and illuminate gas at the nearby contact dis- 
continuity. The IC emission at this time comes mostly from 
high-energy electrons of the FS. The synchrotron emission ex- 
tends to hard X-rays. It comes mostly from the region be- 
tween the FS and the contact discontinuity, somewhat increas- 
ing in intensity towards the CD on account of incre asing MF 
(iRosenau & Frankenfhall[T976HLvutikov & Pohlll2004l) . The ra- 
dio emission is clearly dominated by the emission from the CD, 
where both MF and low-energy electrons are in excess. We 
point out that the decelerating CD is subject to the Rayleigh- 
Taylor instability, which we cannot treat in our 1-D calcula- 
tions. Therefore, real brightness maps would show considerably 
more patchy structure near the contact discontinuity. The ide- 
alized picture presented here, nonetheless, should trace the gen- 
eral trend of the enhanced emission near the CD of core-collapse 
SNRs. 

The Type-IIP SNR at the same age shows quite different 
spectra and morphology of emission. Hadronic emission comes 
from the entire remnant with a significant contribution from the 
CD. The PD spectrum is therefore very soft as the emission re- 
gion is dominated by the low-energy CR population accelerated 
at the RS. It is worth noting that the turnover in the PD spectrum 
is around 10 GeV, a featu re observed in other core-collapse rem - 
nants as well (e.g. Cas A lAcciari et al.ll20Tot lAbdo et al ] l201d> . 
The IC emission is produced in the region between the CD and 
the FS, but its intensity is low. The whole SNR is very bright in 
synchrotron X-rays as the MF is high throughout. This makes it 
possible to clearly distinguish the RS, the CD, and the FS. The 
radio emission, as in Type-Ic case, is dominated by the CD re- 
gion. 

At the age of 1000 years, the gamma-ray emission from 
Type-Ic SNR is dominated by IC radiation. The spectrum of the 
emission is double-peaked: at low energies it is dominated by 
the particles recently accelerated in the shock collision whereas 
the high-energy IC emission is produced by particles accelerated 
at the FS. The PD spectrum is of low intensity with a signifi- 
cant bump at low energies. The IC and PD morphology maps 
at 1 TeV are dominated by high-energy CRs, and therefore no 
hint of emission from RS-accelerated particles is seen. The PD 
emission of FS-accelerated protons illuminating the dense ejecta 
is still visible. The CD is now much brighter than the ejecta re- 
gion, because the target density at the CD is of the same order 
as in the ejecta, but the CR density is much higher. The region 
between the CD and the FS, where the CR density is high, is 
bright in IC emission. Since the MF in the shocked region be- 
came approximately uniform, X-ray synchrotron emission pre- 
dominantly comes from the same region, thus repeating the IC 
morphology. Likewise, the synchrotron spectrum largly resem- 
bles that of IC emission, except for minor differences on account 
of MF radial profile. The radio emission predominantly arises at 
the CD, where the MF is higher than at the RS itself. 

The very high gas density and fairly uniform distribution 
of CRs in 1000-year old Type-IIP SNR renders the PD process 
dominant in the emission spectrum. However, the emissivity is 
mostly localized in the ejecta region. The CD is very bright, and 
the abundance of low-energy protons there makes the PD spec- 
trum soft. The IC emission originates mostly in the weakly mag- 
netized ejecta and the region between the two contact disconti- 
nuities (one behind the reverse shock and the other behind in- 
termediate shock). Both CDs are very prominent is synchrotron 
X-rays. The radio synchrotron map is again dominated by the 
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Fig. 9. Time evolution of intensity maps of Type-Ic (top) and Type-IIP (bottom) SNRs at 1 TeV due to pion-decay (PD) and inverse 
Compton (IC), at 3 keV (SX) and at 1 .4GHz (SR) due to synchrotron radiation. The left column is for the age of 400 years, the 
middle column at 1000 years, and the right column at 1975 years. The scale is linear from zero to maximum in each image. All 
images are normalized by the FS radius, Rf$, at respective time. 



CD behind the reverse shock on account of the distribution of 
low-energy electrons. 

At 1975 years of age, the high-energy radiation from Type- 
Ic SNRs is strongly dominated by IC scattering, with a signif- 
icant contribution coming from low-energy electrons at the RS 
and ejecta regions. Most high-energy IC emission comes from 
the shocked region, where the electron number density is high. 
The PD radiation is very dim and arises mostly in the ejecta, 
with some contribution from the shocked region and the CD. 
Synchrotron X-rays are emitted in the region between the FS and 
the CD, where both the number density of high-energy electrons 
and the MF are high. The radio emission is very bright around 
the CD and the RS, where the density peak of low-energy elec- 
trons is located and the MF is high. Since the ejecta has very low 
MF, the low-energy peak in the synchrotron spectra has nearly 
disappeared. This demonstrates how the radial profile of the MF 
can affect the total emission spectrum, and that a simple scaling 
of the synchrotron emission to IC does not work. 

The brightness of a 1975-year old Type-IIP SNR is still 
dominated by the pion-decay emission coming from the com- 
pact ejecta region and the bright CD, but only up to energies of 
around 1 TeV. At higher energies IC emission takes over, which 
arises throughout the remnant out to the intermediate shock, pre- 
dominantly tracing regions of low MF. In the region of large 
MF, between the intermediate and forward shock, we observe 
hard X-ray synchrotron emission. Radio emission is significant 



in the ejecta region only. One can see that not only the spectra 
of IC and synchrotron emission do not resemble each other, but 
also their morphology is strikingly different. Thus care must be 
exercized in the interpretation of observations if one attempts a 
reconstruction of the high-energy morphology from radio and 
X-ray synchrotron maps of core-collapse SNRs dPetruk et al.l 
2009). A clear shell-type SNR in X-rays may appear as a com- 
pact structure in high-energy gamma-rays, which may be mis- 
taken for a pulsar-wind nebula, thus confusing the identification 
of high-energy gamma-ray sources. Note that the very dense CD 
and ejecta have not expanded significantly, and the radius of the 
CD is only around 3 pc after 1975 years. Resolving this region 
is a challenge for the current TeV-band telescopes, because the 
angular diameter would be only 0.17° at a distance of 2 kpc. 

Alfvenic drift softens Type-Ic SNR spectra rather signifi- 
cantly and renders some of the interesting spectral features less 
vivid, especially in leptonic emission. For Type-IIP spectra the 
inclusion of Alfvenic drift does not have any significant impact 
as discussed in subsection l5.ll 

6. Conclusions 

Using realistic models for the circumstellar medium created by 
stellar mass-loss around massive stars, we have studied the evo- 
lution of the resulting SNR shock waves within this medium, 
and computed the very high-energy signatures from the same. 



11 



I. Telezhinsky et al.: Acceleration of CRs by young core-collapse SNRs. 



Our calculations of the medium into which SNRs evolve take 
into account the changes in mass-loss rate as the progenitor star 
evolved through various stages. In this way, our calculations far 
more precisely identify the various structures in the surrounding 
mediu m as compared to previous authors. For instance, Caprioli 
d201lh assumed a single model to describe the circumstellar 
medium around all massive stars, with a RSG wind region ex- 
tending into a constant density W-R region. Not only is such a 
single description unsuitable over the large parameter space of 
stellar masses, but the structure they have described is not con- 
sistent with observations for massive stars of any given initial 
mass. 

iPtuskin et al.l J2010) also attempted to model the cosmic -ray 
acceleration for 4 types of SNRs, including 3 core-collapse types 
and Type la's. However, their environments were not as detailed 
as ours, and it does not appear that they incorporated the transi- 
tions in density that characterize the circumstellar medium in our 
work. As we have shown, the effects of these transitions makes 
a large difference to the particle spectra and the radiation signa- 
tures. 

We have advance d our model of particle acceleration and 
propagation in SNRs dTelezhinskv et al.ll20 12a b) by adding re- 
alistic transport of magnetic field. Our calculations suggest that 
the magnetic field at the reverse shock is sufficient to acceler- 
ate particles to TeV energies. Although the maximum energy of 
cosmic -ray at the reverse shock is lower than at the FS, their 
intensity is high on account of the high ejecta density in core- 
collapse SNRs. Thus the reverse shock contributes a significant 
fraction of the total CR intensity in SNRs at early times. The RS- 
accelerated particles soften the total emission spectra, in agree- 
ment with recent data. 

Emission spectra from core-collapse SNRs are complex and 
reflect both the distribution of gas or magnetic field and the 
spectral differences between the particle populations coming 
from the forward and reverse shock. High-energy emission from 
Type-Ic SNR is dominated by leptonic processes, particularly 
so when the SNR is older. High-energy radiation from Type-IIP 
SNRs is very soft and strongly dominated by pion-decay emis- 
sion, except for emission at the highest energies from old SNR. 

Both types of remnants show, in addition to a shell structure, 
a center-filled morphology of high-energy emission, because the 
dense and weakly-magnetized ejecta in the interior permits high- 
energy particles to propagate and radiate there. The shell-type 
morphology of Type-IIP SNR, however, may be difficult to re- 
solve with the current generation of gamma-ray telescopes. We 
also find that the complex radial distribution of magnetic field 
causes significant differences in the morphology of synchrotron 
and IC emission from core-collapse SNRs, and therefore a "re- 
mappings" of the spectra and morphology of the two types of 
leptonic emission is not advised. 

The high-energy emission from Type-Ic SNRs is roughly an 
order of magnitude fainter than that of Type-IIP SNRs at early 
ages, but Type-Ic SNRs are good targets for the current gen- 
eration of telescopes at the age of a thousand years and later, 
whereas Type-IIP SNRs lose brightness very rapidly. 

Our work demonstrates that the hydrodynamics of SNRs and 
their interaction with the environment play a significant role in 
shaping particle and emission spectra as well as intensity maps. 
The emission from core-collapse SNRs is significantly different 
from that of Type-la SN R evolving in a uniform environment 
dTelezhinskv et alJl2012ah . 



Appendix 

To derive Eq. |4] we rewrite Eq. [2] in spherical coordinates and 
use the properties of the spherically symmetric flow profile with 
velocity, u: u r - u r , ug — 0, = and du r /d0 = 0, du r /d<f> = 0. 
The magnetic field, B, has all three components, B r , Bg, B$, and 
V • B = 0, so 

u(V • B) = (9) 



B(v ■ u) = 

B(^ + 7ibl (*«M> + 7^|f)= (10) 
B r (¥ + fr) r + Bg (2* + + + £ 

(B • v)u = 

[ R ihk. _j_ El _i_ Blt> dj*r _ B$uq+B^U£ \ a 

\" r dr + r 36 + rsin(fl) 8$ r )' + 

in due . B e dug , B <t Sue , B g u r _ Bgi^cot(fl) \ % 

\ r dr r d6 r sin(0) dip r r ) v / 

In d"f , Bu £«£ B# ditj B#u r B t u„ cot(fl) ^ 1 _ 

\ D r dr + r 86 + rsin(fi) d<t> + r + r ) 9 ~ 

B r %t + Bg±d + B#i 



(u • v)B = 

/ dB r uedB, . "4, dB r _ ugBe+u^ \ « 

\ Ur dr + r aB + rsin(fl) 30 r ) ' + 

I dB„ ue_dBe «g BB U u u B r cot(fl) \ g 

\ r dr r d6 rsin(8) dtp r r J \ ' A* 

/ dB^ u^dB^ uj, 8B t u^B,. u$B$ cot(fl) \ " _ 
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Combining the expressions we find for the evolution of the 
three components of B 

dB r D Uj_ n du r . d du r _ dB r 

IT ~ ~ LCS r— ~ D r -^r + Br-fr ~ U r -fr 

= - fr (B,u r ) - 2B r l jr + B r "a7 . 

dB S _ O D "r D dUr , r> U r ,. dB S _ d / D „ \ D «, 

— - ~2Bg- - Bg-^r + Bg- - U r ^r - ~ Tr (BgU r ) - Bg- , 

— OR", d du r D «, dB$ _ d In ,.\ n u r 

-gf - -^«07 - «0"57 + «07 - Ur-gp - -g- r {tifiU,.) - tlrb- . 

(13) 
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